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Abstract 

This work proposes a model for granular deformation that predicts the stress and velocity profiles 
in well-developed dense granular flows. Recent models for granular elasticity (Jiang and Liu 2003) 
and rate-sensitive fluid-like flow (Jop et al. 2006) are reformulated and combined into one universal 
elasto-plastic law, capable of predicting flowing regions and stagnant zones simultaneously in any 
arbitrary 3D flow geometry. The unification is performed by justifying and implementing a Kroner- 
Lee decomposition, with care taken to ensure certain continuum physical principles are necessarily 
upheld. The model is then numerically implemented in multiple geometries and results are compared 
to experiments and discrete simulations. 



1. Introduction 



A general constitutive law for granular materials, valid over all common regimes of deformation, 
remains an important open problem for several communities. Mechanicians and materials engineers 
have modeled granular materials for centuries guided by principles of continuum solid mechanics 
63,81, 3^ Zlj 3i 28 , [ill , where various failure and yield postulates have been employed. In recent 
years, a resurgence of interest in the topic has arisen among physicists [4l|, [11, [4^, [2^, H^, [H| , where 
modeling efforts have been directed primarily toward statistical and fiuid dynamical approaches. 

Some of the difficulty in modeling dense granular materials is due to the existence of highly 
disparate responses when material is flowing "slowly" or "moderately" (see section I4.2|) , both of 
which only vaguely relate to the completely static response. To describe certain phenomenological 
behaviors, an astounding array of physical hypotheses have been suggested. These include: rolling 



and static phase interactions 112 



I. fisl [i3l Bagnold rheology generated b y "g ranular eddies" [2£^ 
granular temperature-dependent viscosityfTSj , density-dependent viscosi ty I57l .[l0j , non-local stress 
propagation 13, U | , free- volume diffusion opposing gravity [HI, [H, [13, [1, l7d|. "shear transformation 



zones" coupled to free- volume kinetics [5J, [55 
Landau order parameter Such theories can typicall 

usually for a specific geometry: avalanc hing surface flows (l^.ll 
Couette cells 57, [13], inclined chutes IjOl [Hol. or wide silo [H 
granular flow theories can be found in 



and partial fluidization governed by a Ginzberg- 
fit certain experimental data 1 59] , 
"l^. [ll , inclined planes [1, [23 , 
164. Isl. [t^I . A good review of 




Some models have been developed that a re g eneral enough for multiple flow environments [l|, [t^. 



47], IJ]. In particular, the work of Jop et al. (4j| introduces and tests a generic, 3D, non-Newtonian 



fluid constitutive law representing "moderate" sized granular shear-rate as a local function of the 
Cauchy stress. The material is treated as a Bingham fluid and a computation of the stress and 
flow profiles can be made in regions above a Drucker-Prager yield criterion. While a breakthrough 
by many standards, the law makes rigid-body assumptions in static phases where stresses are not 
computed. In many circumstances, granular assemblies can establish static zones at walls (such as 
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in flowing hoppers and silos [63, [Hi), and effective flow computation is impingent on the ability 



to model traction-based interaction laws between sub-yield material and the system boundaries. 
Stress-based quantities in sub-yield zones are needed to approach a variety of open questions in 
granular flow theory, such as how "kick-stresses" [6^, |2l[ develop in draining storage bins. 

This work seeks a general constitutive law for 3D, dense flowing granular materials. Specifically, 
we focus on flows with significant regions of moderate shearing, as would be expected in day-to- 
day environments like fiowing chutes, hoppers, and heaps. To this end, the model of Jop et al. is 
integrated with an appropriate granular solid response, by way of a joint elasto-plastic treatment. 
The unified formulation can be used to predict the stress and flow proflles everywhere, for an 
arbitrary geometry with arbitrary admissible kinematic/traction boundary conditions. For sufficient 
generality and accuracy in describing the static phase, we ultimately choose to pair with the granular 
elasticity law proposed recently by Jiang and Liu [44] . This happens to be a natural choice for 
several reasons, not the least of which is its proven capability to represent several static and acoustic 
phenomena observed in granular matter. 

The process of constructing a combined law is itself non-trivial. Both components require some 
degree of reinterpretation — the Jop et al. law is converted to an evolution law for the plastic 
deformation gradient, and the Jiang-Liu elasticity law is modified to permit flowing states while 
maintaining elastic deformation. Moreover, some essential postulates of the chosen elasto-plastic 
framework are not obvious in well-developed flow, requiring additional substantiation. 

This paper is organized as follows. First, we describe some general requirements for a continuum 
treatment of granular statics/flow (section [5]). Then a discussion of past solid and flow constitutive 
laws for dense granular media are reviewed (sections[3]and[?]), ultimately culminating in a description 
of the two laws mentioned above. The process of combining these responses is discussed (section 
[5|), which ends with a summary of equations for the proposed model. Finally, section [6] provides 
discussion on the numerical implementation of these equations and displays results in several flow 
geometries, with comparisons against the data of 44|,[5^, and [t^. 



2. Granular matter as a continuum 

Before development ensues, we take a moment to lay out the limits of validity for a first-order 
continuum treatment of a discrete granular collection, as shall be proposed herein. For a more 
detailed set of arguments on deterministic limits in random composite materials, the reader is 
directed to the seminal work of Ostoja-Starzewski [67.] [66] . 

Recent Discrete Element Method (DEM) simulations of [7^ have shown that in multiple 3D well- 
developed fiow environments, a dense granular element of width 5d (for d — particle diameter) 
captures many of the plastic flow properties expected of an RVE. Among neighboring volume 
elements of this width, the average stress and deformation rate appear to vary smoothly. Within a 
fiowing element, the eigenvectors of the instantaneous space- average Cauchy stress align to a high 
degree with those of the deformation rate tensor, evidencing the onset of a deterministic relationship 
between the two fields. Moreover, a predictable dependence of the packing fraction on the pressure 
and shearing rate emerges at this length-scale. Of course, the representative element width is not a 
universal constant; environments with large spatial gradients necessitate smaller elements. Thus the 
5d element width can be interpreted as a common limit for day-to-day fiows, compacted primarily 
by gravity, where boundary conditions on the global body vary on a scale larger than O{10d). 

The result is not outlandish in light of the past observations by Glasser and Goldhirsch [2^ 
and Goldenberg et al. [sij, where in depth 2D studies on disk arrays were performed to quantify 
the effects of spatial averaging on granular stresses. Rycroft's result is in line with the extremal 
trends observed over the spectrum from rapid/dilute disk flows to static/elastic deformation of a 
disk assembly. Furthermore, it is found that static materials gain resolution independence on a 
similar scale to flowing materials, suggesting that a 5d element width could give valid RVE behavior 
for dense granular materials in either static or flowing states. 
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Previous arguments against a granular continuum, particularly in the static phase but also 
during flow, have focused on the presence of force chains in granular matter. In 2D experiments 
and simulations of disk assemblies |2^, concentrated chains of interparticle force have been 
shown to exist over many particle widths. It was argued that if forces are not homogeneous at a 
meso-scale, continuum relations at this scale cannot exist. However, in 3D simulations of flowing 
monodisperse sphere packings, we observe that the force chains have a dramatically shorter length 
(76| . A possible geometric explanation for this phenomenon is that a 3D granular assembly has a 
much higher average coordination number, reducing the likelihood that only two contacts maintain 
the majority of the force on one grain. These simulations also include interparticle contact friction, 
which may contribute to the smoothing out of force chains as has been previously shown in 2D 



static disk assemblies [34|. Typically speaking, whether flowing or static, we observe that a 5d 



granular volume element contains a diffuse network of contact forces, enabling a sufficient degree 
homogenization of the stresses at this scale. 

Frequently, a granular collection that has been set into motion will possess regions of fluid-like 
flow adjacent to essentially solid-like regions. Examples of such flows include wide draining silos 
[lil . [t^I and hoppers (so-called "core flow" ) where a broad column of material extending upward 
from the orifice flows like a fluid, while regions closer to the side walls remain almost completely 
static. As will be discussed later, perfectly clear solid/fluid interfaces are rarely observed in granular 
flow, which has led some to believe that the solid-like zone is actually just a 'highly viscous' fluid 
region [57I [lo| . Solid- like material does undergo intermittent rearrangement events when close 



to a zone of moderate flow-rate, but we flnd that the stresses in these regions have essentially zero 
rate-dependence. For example, when a DEM simulation of silo flow is suddenly arrested, say by 
shutting the oriflce, the stresses in the solid-like regions remain virtually unchanged, supporting 
static shear stress like a solid. It appears that a mechanism disconnected from the flow-rate or any 
notion of viscosity is responsible for maintaining the stress tensor in solid- like granular matter, even 
if occasional failure events are occurring within. This suggests that a continuum description for 
granular statics is likely to be disjointed from the flow description. Hence, we seek a general model 
that splices separate statics and flow laws into one constitutive form. 



3. Continuum statics 

3.1. Stress-only laws vs. elasticity 

To close the equations for static force balance, two approaches have been utilized: stress-only 
laws and elasticity. Stress-only relationships constrain the stresses directly, by asserting that the 
stress tensor must satisfy some a priori relationship. Examples include: Janssen's law of differen- 
tial slices (originally proposed by Janssen in 1895) where vertical and horizontal stresses are set to 
be proportional, limit-state Mohr-Coulomb plasticity or the "Incipient Failure Everywhere" (IFE) 
hypothesis [si] where a Mohr-Coulomb failure line is required to exist at all locations, and "Ori- 
ented Stress Linearity" [ul where stresses propagate in directions aligned with the presumptive 
microstructure of the packing. 

While stress-only relationships are convenient and have had some success, their physical assump- 
tions can be questionable. For example, static granular matter is rarely in a limit-state of incipient 
failure [76] and wall shear is not compatible with a Janssen-style stress tensor [63] • Most stress-only 
laws are defined only for 2D media, which brings out issues of generality. They often predict a 
"hyperbolic" character to the stress profile, where stress quantities propagate in certain directions 
from the boundaries. The notion of force propagation was backed chiefly by the observation of a 
double-peak pressure distribution beneath a bed of grains undergoing a point force from the top. 
Work by Goldenberg and Goldhirsch 32] has shown, however, that in the presence of interparticle 
friction and a large system size to particle size ratio (as commonly found in engineering applications) 
the pressure distribution is indeed a broad single peak beneath the point force, as one would expect 
for an elastic bulk media. 
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This brings us to elasticity, which shall be our preferred method for granular statics. It is a 
sensible approach seeing as the grains themselves are elastic bodies presumably enabling general- 
ization to an Effective Medium Theory (EMT) where grain-level elasticity extends statistically to a 
continuum mean- field theory 22, 2^, S^, H^- Reversible (elastic) deformations have been observed 
in granular matter for strains less than 10~^ [5l[. This is negligibly small compared to the size of 
typical plastic deformation. However, grains are commonly composed of stiff material indicating the 
important role that small elastic strains may play in the development of the stress profile. It could 
be argued that the barely noticeable elastic strain of a static assembly is what originally impelled 
scientists to seek stress-only laws, so that rigid-body assumptions could be used. Rather, it seems 
there is no generally applicable stress-only constraint to accurately define a 3D stress tensor. The 
elastic strains are small but non-ignorable, and bear essential importance to the stress description. 

Bulk elasticity for cohesionless grains is not likely to have a simple form since, for example, the 
material is unable to support tension and thus the small strain response cannot be approximated as 
linear. Nonlinear EMT has been derived from Hertz-Mindlin interparticle contact mechanics 22, gH 
and modified by others [i^, offering reasonable albeit not completely satisfactory agreement with 



experiments [51l |30 | . A recently proposed elasticity model encompassing many of the same features 
as EMT was proposed by Jiang and Liu [i^l in 2003. This state-of-the-art formulation has had 
multiple successes and is well-suited to our end goal of combining with a plasticity model. 



3.2. Effective Medium Theory for hulk granular elasticity 

In the classical work of Hertz, two perfectly elastic spheres that are pressed into contact repel 
each other with a contact force that depends on the radii of the spheres and the "apparent overlap" . 
For two spheres of identical radius R located at xi and X2, the normal force contact law is 

where 5 is the normal overlap (1/2) (2i? — ||x2 — xi||). The parameter fc„ is an effective stiffness 
equal to 4Gg/(l — Vg) where Gg and Vg arc the shear modulus and Poisson ration of the sphere 
material. 

Presume a granular collection with an isotropic distribution of contacts that undergoes affine 
deformation when stressed. The Hertzian contact law generalizes to an EMT that gives the following 
bulk modulus under pure compression: 

At = — $Zfc„(-itrE)i/2 

where $ is the volume fraction, Z is the mean coordination number, and E is the (infinitesimal) 
strain tensor. The nonlinearity can be seen as arising directly from the nonlinearity of Hertz's 
contact law. The major point is that the bulk modulus scales with (compressive strain) or 
equivalently as p^^^ under isotropic compression. This has been verified directly in large-scale DEM 
simulations of compressed sphere packings both with and without interparticle contact friction jsl] . 

Supposing frictionless spheres, a mean- field shear modulus can also be derived in a similar fashion 
by analyzing an arbitrary affine deformation including shear strain. One ultimately finds that the 
bulk and shear moduli under non-zero shear scale similarly to the above form: 

K cx G cx $Z(-trE)^/2 ^ (^zf/^p^^^ (2) 

The inclusion of shear strain renders EMT less accurate. The moduli agree with the above scal- 
ings but only under low pressures (<10 MPa) Experiments verify that the shear dependence on 
p is characterized by a larger exponent at higher pressures [30| . This could be because shearing un- 
der higher pressures tends to render the affine displacement assumption less valid. Time-dependent 
relaxation occurs, which significantly complicates a determination of the shear modulus. Attempts 
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to improve the EMT elastic moduli during shear are not vastly affected by including tangential 



forces in the analysis; the inclusion of Mindlin's tangential force law [6l| between spheres does not 
change the above scalings. 

3.3. The Jiang-Liu granular elasticity law 

In 2003, Jiang and Liu proposed a modified approach to granular elasticity Rather than 
continue laboring on a microscopically derivable mean field theory, their aim was to augment the 
more successful results of EMT with presumptive forms that capture known macroscopic behavior. 

They consider an elastic free energy density of the following form: 

V'(E) = SxM:QA2+7Ve) (3) 

where S is a relative stiffness that can vary with packing fraction. The compressive strain and shear 
strain are measured respectively by 

A = -trE, 7=|Eo| 

where for any tensor A, we adopt the notation Aq = A — (l/3)trA and |A| = y'X^ij- ^fj- The 
dimensionless parameter ^ will be discussed shortly. In the small displacement limit this gives the 
Cauchy stress 

G(A) k(A,7) 

Under isotropic compression, the pressure is proportional to A'^/^, in agreement with the most 
successful result of EMT. The shear modulus scales with IS^/'^ as in EMT, but the added nonlinearity 
in the full form of the bulk modulus allows for some important properties not well-represented by 
EMT. 

In an arbitrary state of strain (7, A), the corresponding Druckcr-Pragcr friction state is 

^ = !£i = V2 ((7/A)-ie + (7/A)0.5)"' 
p 

By varying the input r = 7/ A, one can show using basic calculus that /i has a global maximum 
value of Thus, by selecting f accordingly, we can prevent certain states of friction from ever 

arising. This is one of the defining features of Jiang-Liu elasticity. 

In the work of Jiang and Liu, static materials were the primary concern. Hence, a static yield 
criterion was declared and ^ was set to /ij^ so as to require that no purely elastic state can exist 
above /i^. The work herein shall attempt to integrate elasticity within a complete elasto-plastic 
framework, so the basis for selecting ^ will be slightly different. At the point where the stress 
state reaches fJ-max, the elasticity relation loses convexity. The particular condition for convexity is 
|Eo| < (p/2_B)^^^ a/2C- Figure [T] elucidates the fashion by which increases to its maximum under 
increasing shear strain. Plots are discontinued in non-convex regions, where the elasticity law is 
deemed void. 

More than just a model that connects with yield, the Jiang-Liu elasticity model has had some 
convincing successes in representing granular statics: 

1. The nonlinear form of T produces a stiffness tensor g^j^ that agrees to a satisfactory extent 
with the form of the stiffness tensor extracted from experimental data (5l] |. 

2. The model predicts Janssen-type saturation of the wall stresses in a tall silo. The ratio of vertical 
to horizontal compressive stresses in the silo is found to be approximately constant when not 
close to the walls. This verifies the commonly used notion of a "coefficient of redirection", 
which has been verified in DEM simulations [76|. 



5 




0.4 0.6 0.8 

v(2B/p)=^(29-"^ 



Figure 1: The shear stress vs. shear strain relation under Jiang-Liu elasticity represented as a single dimensionless 
plot (left), and plotted equivalently as a family of shear stress vs. shear strain curves in SI units using ^ = 4, 
B = 7 X 10^ Pa, each curve determined by the applied compressive pressure (right). 



3. The model predicts that a granular material under simple shear stress responds anisotropically 
to the addition of a point-load at a surface. Such anisotropy under pre-stress is a well-known 
granular phenomenon that is captured appropriately in the nonlinearity of the Jiang-Liu model. 

4. It has been observed that preparation history is largely responsible for the stress dip that one 
often observes under the peak of a sand pile. This fact is reproduced by the elasticity model 
when solved assuming an initial packing fraction distribution that one might expect for a conical 
pile constructed from grains flowing out of a nozzle. 



4. Continuum flow 

4.I. Bagnold scaling and relevant dimensionless quantities 

We now move on to a discussion of granular flow laws for inelastic deformation. Bagnold's 
seminal work on granular flow followed from shear experiments on granular/fluid suspensions. In 
the "grain-inertia" regime where the effects of the interstitial fluid are small, Bagnold found that 
both shear and normal stresses on the shearing wall depend quadratically on the shear rate [6J, a 
phenomenon that came to be known as "Bagnold scaling" . Bagnold scaling has been verified for 
dry grains in a number of experiments and simulations [56|, l8(i [ssl [tiI [17| . An explicit form in the 
case of simple shear of a dense configuration of dry grains can be expressed as the following pair of 
dimensionless relations: 



$ = for/=^^ (5) 

M = 9(1) for M = ^ (6) 

In the above, $ is the packing fraction, P is the pressure on the shearing plate, r is the shear stress, 
and the steady shear rate is 7. The dimensionless number / is commonly referred to as the inertial 
number or normalized fiow rate, and fi is the effective friction. 

The simplest way to understand these equations is through dimensional analysis. The major 
physical quantities involved in a gravity-free simple shearing between long rough plates are the 
material parameters d and ps, and the variable quantities $, t, P, and 7. This tacitly ignores the 
possibility of any other length-scales playing a role and presumes that collisions are fully dissipative 
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Figure 2: Qualitative diagrams (primarily due to 17]) showing the variation of dimensionless parameters through the 
various flow regimes under simple shearing. In other geometries, the quasi-static regime of fi is not as clear to define 
in terms of I due to meso-scale effects. The coefficient of restitution affects <& vs. / in the coUisional regime. The 
moderate regime is relatively well-determined as a sole function of /. 



(pressure high enough to damp out restitution), two assumptions whose consequences are important 
and wiU be discussed shortly. The particle-on-particle contact friction /ip is also ignored. Granted, 
fip does affect fi, but it has been found to merely translate the fj, vs. / relationship vertically [it} . 

With these assumptions, the packing fraction and effective friction in a simple shear apparatus 
should arise uniquely from P and 7 as a matter of cause and effect. Nondimensionalization implies 
that fj, and $ should depend only on /. Unlike a Newtonian fluid where a temperature time-scale 
exists, the quadratic dependence of stress on shear rate can be seen to arise from the fact that 7 
can only be scaled by the square root of a stress quantity. 

4-2. Flow regimes 

It is important to clarify the properties and overall validity of equations [5] and [S] over the range of 
possible flow behaviors. First, zero coUisional restitution was presumed in justifying these equations. 
The assumption is valid if the kinetic energy of a collision is always dissipated in full on impact, 
presumably in the form of heat and sound. 

For faster shear rates, this source of energy loss causes notable rate-sensitivity. Resultantly, / 
becomes one-to-one with /i, and 7 is immediately determined by r and P as in a non-Newtonian 
fluid. Gathering these properties into a general classification scheme, a flow rate is deemed to be 
moderate when / is large enough for rate-dependence, but small enough for the flow to stay dense 
as per the coUisional collapse argument. Data of da Cruz et al. would suggest this regime 
lies within the band 10~^ < I < 10~^. In day-to-day terms, the flowing region of an hourglass is 
typically in the moderate range. 

Moderate flows have the property of shearing dilation, where increasing the normalized flow rate 
causes the steady-state packing fraction to decrease (i.e. / in equation [5] becomes a decreasing 
function). This should not be confused with shear dilation, which refers to the drop in packing 
fraction as a function of total shear that occurs to a dense assembly at the beginning stages of a 
shear deformation. Flows too slow to be deemed moderate may still undergo shear dilation due to 
packing geometry, but rate effects like shearing dilation only set in for faster flows. 

Beyond the moderate range, dilute or coUisional flows occur in general for / > 10^^ and cor- 
respond to the breakdown of the zero restitution assumption. When / becomes this large, particle 
collisions are accompanied by some additional "bounce-back" akin to a gas. The collisions are chiefly 
binary, and particles rarely maintain long lasting contacts. These flows may require a temperature 
quantity to store fluctuational energy and can be described with dissipative Boltzmann kinetics. A 
sand blast could be considered a basic example of the coUisional regime. 
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On the other side of the spectrum, where / < 10~^, we enter the quasi-static regime. The pack- 
ing fraction in simple shear does not vary noticeably with / in the (time-averaged) steady limit — 
the inertial time is always small enough for the particles to find tight compaction. Without a sig- 
nificant contribution from coUisional dissipation, rate-dependence subsides, and more complicated 
dissipation mechanisms dominate such as frictional contact sliding. The stress/strain-rate relation- 
ship becomes singular; driving the system with a range of quasi-static normalized shear rates all 
give the same time-average value for fj,. Soil creep is typically in the quasi-static range. 

The discussion thus far has focused on simple shear, where stresses and time-average flow are 
spatially uniform. While the other regimes display a strong local rheology, quasi-static flows are 
sensitive to grad ients in the flelds as expressed through some non-local term with an associated 



length-scale [59|, |8|, |47|, Ij. For example, consider steady flow in an annular Couette cell. Slow 
moving material far from the inner wall is observed to constantly creep [H^, even though the 
stress state should be below yield by all common local measures. Indeed, the motion appears to 
be caused by a non-local effect where faster flow near the inner wall has effectively "bled out" 
into neighboring material. Grain-level specifics such as roughness, grain shape, and conflgurational 
statistics (including wall effects) should affect the non-local flow behavior through the new length- 
scale. The size, dynamics, and general interpretation of the length-scale are object of debate, though 
most agree its size should be on the order of several particle widths. 

Observe figure[5]for a schematic view of how /i and <I> vary throughout the flow regimes (in simple 
shear). The major points of this discussion are: moderate flow is much simpler than quasi-static, 
and coUisional flow is outside our interest as it is not dense. Moderate flows are characterized by a 
local rheology relating / to both fi and <&. The regime is definitively rate-dependent and equation [5] 
inverts into a fluid-like law wherein the flow rate can be determined uniquely from the stress state. 
As / decreases to the quasi-static limit, however, quantities that were previously ignorable become 
important. Grain motion is locally correlated at some length-scale causing nonlocality and a larger 
role of grain-level properties. The flow does not permit a fluid-like treatment as before, since the 
dissipation is largely rate-independent. 

4-. 3. Quasi- static flow models 

To review dense flow models, we begin with those aimed at quasi-static behavior. Since the 
stresses are not a direct function of the flow rate, they are typically obtained by asserting an 
elasticity law upon a small component of the total deformation. Linear elasticity is often presumed, 
but this may be an oversimplification that carries more consequences at lower stresses as described 
in section[3l Critical State Theory f79| . R udnicki-Rice-tvoe modeling the Anand-Gu model [H, 
and the model of La Ragione et al. |72l | each can be used for granular deformation in this fashion. 
Other models couple to the IFE stress formulation in 2D, such as rate- independent coaxiality [6^ 
and the Stochastic Flow Rule . 

Some quasi-static theories account for nonlocality, though none have been proven to do so reliably 
for steady flow in arbitrary 3D environments. Several theories are based on new definitions of tem- 
perature, which introduce a length-scale via a "heat equation" . These include Shear Transformation 
Zone (STZ) Theory 0,1231' which relies on an effective temperature governing STZ creation, the 
dense fiow theory of S. B. Savage (tsI, w hich defines a granular temperature to measure strain-rate 
fiuctuations, and Edwards statistics 23], which utilizes a temperature-like 'compactivity' derived 
from an entropy per free volume. These models provide interesting physical insight, but do not 
appear to be at the point of development that simulating arbitrary fiows would be possible — some 
are restricted to 2D, the boundary conditions for the new temperature are rarely obvious, and the 
equations may not be closed except in a few symmetric test cases. Besides temperature approaches, 
nonlocal behavior could also be described with a diffusing order parameter, as in [iJij or through 



a more general strain-gradient plasticity theory [33|, |3J, |88|, |38|, |37 

The quasi-static flow regime, though important, appears at the moment too difficult to account 
for appropriately within a simple 3D continuum framework. Likewise, the following concession 
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Figure 3: The plastic flow rheology in simple shear plotted as one dimensionless relationship (left), and plotted 
equivalently as a family of shear-rate vs. shear stress curves in SI units, each curve determined by the applied 
compressive pressure (right). 



is enacted: The model to be constructed in this work shall neglect quasi-static flow 
behavior altogether, opting instead to combine a static response directly with moderate 
flow rheology. The model herein may ultimately serve as the backbone for a fuller model that 
also incorporates the dependence of a length-scale on the slow dynamics. This possibility shall be 
considered in more depth when we compare predictions of the model directly to experimental/DEM 
data. But for now, we accept inaccuracy in describing quasi-static motion in exchange for a closed, 
general model capable of giving worthwhile predictions over the full range of dense material behavior, 
accounting for both statics and flow. 

4- 4- Moderate flow law of J op et al. 

A closed form law to predict moderate flow must now be discussed. Since the moderate flow 
regime is (monotonically) rate-dependent, the function g in equation [S] should be invertible. Hence, 
increasing the normalized shear rate / results in higher /i. This notion may seem counterintuitive 
to the fundamental idea that ^ should decrease to a kinetic value as the rate of sliding picks up. 
But recall, for moderate flow rates, the impact dissipation dominates sliding affects. In slower or 
transient flows, shear weakening is indeed observed for "overconsolidated" material, and accounted 
for in various models via hardening parameters [J, . 

Based on results from numerical simulations of planar shea r [TtI . |40| , and successful extensions 
to plane-strain inclined chute flows 5^, S^l, the experiments of [43j were conducted to quantify 
for glass beads: 

/ = g-i(/i)=/o^^^ forM>/i,. (7) 

The values of the parameters were measured at /q = 0.279, /i^ = tan 20.9°, and /i2 = tan32.76°. 
The relation states that the normalized shear rate / increases as the material is sheared with higher 
^. But /i must exceed some static yield value before any plastic flow ensues. There is also some 
maximal /i value, 1x2 ^ and all steady shear flowes should be tenable for a value of /i less than /i2. 
Consequently, if applied stresses exceed /i2 the flow is deemed accelerative with no apparent steady 
state. While this last point may be somewhat of an oversimplification from a physical point of view, 
the approximation is helpful especially in light of the fact that /i for dense flow is usually well below 



The first attempt to extend equation [7| into 3D was met with high success. In Jop et al. j44l |. 
codirectionality was applied, presuming that the deformation rate tensor 

D = (1/2) (Vv + (Vv)^) 
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is proportional to the dcviatoric stress tensor 



To = T-(l/3)(trT)l. 



Written in full, Jop proposed the following generalization of equation [71 



D 



d 



fj. - /is To 



(8) 
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where now P — — (l/3)trT and ^ ~ t/P where r = |To|/-\/2 is the equivalent shear stress. When 
/i < /is , we take D = establishing a Drucker-Prager yield criterion. Since the flow condition being 
used is codirectionality, the flow rule is thus non-associative. 

Since D is proportional to a deviatoric tensor, the flow rule asserts plastic incompressibility. 
While, as previously noted, dilation in dense flow does occur, it is typically on the order of only a 
few percent and quickly reaches a steady value over large deformations. Moreover, <I> is unnecessary 
to compute the steady shear rate, hence the approximation of plastic incompressibility should have 
negligible effect on the velocity field of a dense steady flow. Note that evolution laws for the packing 
fraction have not yet been quantified in this context. Some quasi-static flow models attempt this 
Qjllijlzll but rather than try to modify one of these, we go along with the presumptions of Jop et al. 
and ignore plastic dilation. The bigger impact of this assertion is not on the flow, but rather on the 
stresses in the static regions — If plastic flow transiently passes through a region that is ultimately 
static in the steady-state, the dilation that occurred there affects the local elastic moduli. 

The flow model was tested experimentally in a chute apparatus with rough sidewalls and bottom 
to induce non-trivial 3D flow. Results were shown to match experiments to a high degree (always 
within 15%) even while varying several parameters (e.g. inclination angle, flow height and width). 
With positive experimental validation, the flow model became one of the first to describe well- 
developed, inhomogeneous, 3D, dense granular flows under moderate flow rates. 

We now proceed to unite the Jop et al. plasticity law (equation [5]) with the Jiang-Liu elasticity 
law (equation |4]). 

5. Combining elasticity and plasticity 

5.1. Kinematics 

First, some kinematic quantities of interest to the discussion shall be briefly reviewed (see [s^ 
or for more details). For a material element that begins at location X and resides at x some 
time t later, the motion Junction x is defined by x = x(X,t). The elastic and plastic responses 
are combined using the Kroner-Lee decomposition of the deformation gradient F = dxO^Tt) / ffX. 



The decomposition models the deformation of an element at any time as a progression of two 
stages of deformation. First, the unstressed reference configuration is deformed by the plastic 
deformation F^ to an intermediate configuration, which is still deemed as stress-free. From there, 
the material deforms elastically via F*^ to its final or deformed configuration. The intermediate 
configuration can be thought of as the residual state that occurs when all stresses on a deformed 
volume element are released, thereby unloading all elastic mechanisms, leaving only the total plastic 
deformation. 

The spatial velocity gradient can be expressed in terms of F via 
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where ' always represents the material time derivative. Generalizing this, the elastic and plastic 
velocity gradients are defined as 

LP = FPFP-^ (11) 

which in turn enables definitions of the elastic and plastic deformation rate and the elastic and 
plastic spin: 

D^^ = sym(L"=) = sym(L^') 

W = skw(L'=) WP ^ skw(L^'). 

The deformation gradient and its elastic/plastic counterparts are presumed to have positive 
determinant, permitting the polar decomposition 

F"^ = R'^U'^ (12) 

where XJ'^ is the elastic right stretch (which is symmetric positive definite) and R*^ is the elastic 
rotation. The plasticity model we wish to use is incompressible and as such we have detF^' = 1. 
Consequently, detF*^ = detF. From the elastic stretch, we invoke the material Hencky strain measure 
to define the elastic strain: 

W = log(U'=) (13) 

5.2. Elasto-plastic constitutive picture 

The constitutive framework to be used herein is analagous to a spring/damper combination in 
series, with Jiang-Liu elasticity representing the spring deformation and the flow law of Jop et 
al. representing the damper (with a yield criterion). Notably, the stresses are modeled as being 
determinable everywhere from the field E"^. When the resulting stresses satisfy yield, the plastic 
flow law determines the flow-rate L^. 

It is a non-trivial physical assertion to model all stresses in a dense flowing granular material as 
being derivable from the elastic deformation of the grains. Other microscopic stress agencies exist 
such as internal viscous damping and fluctuational stresses. To validate our assertion, discrete par- 
ticle simulations of Rycroft |74| were performed in several dense, inhomogeneous flow environments. 
The instantaneous stress tensor over an element is represented by the particulate formula 

\ 

r(y) ^ f (y) _ ^(i) - v) (v(*) - v) (14) 

for N the number of grains in the element, V the element volume, f the contact force of particle i 
on particle j, r*^'-'-' the vector connecting the centroids of particles i and j, m*^*-' the mass of particle i, 
v^*) the velocity of particle i, and v the average velocity over all N particles. The first term inside the 
parenthesis represents a solid-like stress response derived from contact forces, while the latter term 
is gas-like, accounting for stresses due to velocity fluctuations. Rycroft found that the gas-like term 
is always exceedingly small compared to the solid term. The contact law used in the simulation 
was visco-elastic (plus sliding friction) giving the force decomposition f^'^^ = f ^^^^ -|- fv(*i) gy 
comparing the viscous and elastic force contributions, Rycroft found that the elastic contribution is 
always vastly dominant. Altogether, 
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This evidence suggests the elastic mechanism accounts for virtually all stress in a flowing or static 
element of dense granular material, lending reasonable validation to our proposed elasto-plastic 
treatment. 
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5.3. Modifying the Jiang-Liu elasticity law 

The Jiang-Liu elasticity law requires some slight modifications in order to be used for elasto- 
plasticity. Since flowing materials can have above /is, the elastic law must permit these stress 
states. Hence, we propose the following important modification to Jiang-Liu elasticity: Set ^ = ^12^ 
instead of ^ = /^J^. Little loss in representing elasticity occurs from this change. The ^ parameter 
was engendered to represent the macroscopic repose angle of a static granular assembly; it was not 
determined from any quantitative microscopic requirements. 

With ^ = /i^^, there can never be an elastic stress state that has fj, > fj,2- Looking back to 
equation[8l notice that the flow rule cannot apply to a stress state that exceeds /i2. Since the Jiang- 
Liu law admits a cap on the value of /i, this property prevents elastic stress states from entering 
the forbidden regime of the Jop et al. flow law, which offers several benefits when attempting to 
numerically solve for flow. Analytically speaking, capping the elastic stresses at /i2 is not necessary 
for a solution; a cap above /i2 would also be acceptable, though the plastic response would always 
preclude such states from arising. Regardless, ^2 is a natural choice that minimizes the extent of 
alteration of Jiang and Liu's original formulation. 

5.4. Mathematical specifics 

With the motivation provided, we now go about providing a mathematically rigorous framework 
to unite the proposed elastic and plastic responses. The following is based on a form for thermody- 
namically compatible elasto-plasticity developed in . An abbreviated discussion shall be provided 
here. More details can be found in [46| . 

At the outset, Newton's equations of motion must be upheld for force and torque balance, and 
mass is conserved: 

^•T + pg = pv, T = T^, p = po(detF)-i (^5) 
ax 

where po is the initial material density, equal to random-close-packing (63%) times ps. To institute 
simultaneous elastic and plastic constitutive responses, three major restrictions are enforced 

1. Frame-indifference 

2. Non- violation of the second law of thermodynamics 

3. Isotropy of the reference and intermediate configurations. 

The presumption of isotropy, particularly in the intermediate state, deserves some attention 
Anisotropy of plastically deformed granular material has been studied extensively both in 2D 



and in 3D where it has been quantified [87| and modeled JM), 86] often by inclusion of a fabric tensor. 
These studies show varying levels of anisotropy that evolve over a range of strains that are large, 
but always significantly less then 50%. In the limit of well-developed flows, however, experimental 
evidence of [SSJ and discrete simulations of 1^, 7^ support the use of an isotropic flow law. 



With the presumed F'^F^' decomposition, one can show via power conjugacy that the stress 
interacts with the elastic and plastic mechanisms through a surrogate known as the Mandel stress 

M = J^F^^TF^ (16) 

To connect to more common stress measures, M is equivalent to C*^Tii where Tn = J'^F'^^^TF'^^-'" 
is the second Piola stress as measured using F*^ instead of F, and = F^'^F'^ is the elastic right 
Cauchy-Green tensor. As in the analogy with a spring and damper in series, the Mandel stress 
determines both the elastic strain and the rate of plastic flow. Equations |3] and [5] are written in this 
parlance to give the elasto-plastic constitutive relations: 



Elasticity relation: 

M = 2GEJ5 + K (trE'^) 1 (17) 
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where k = BVA[1 + -f^/{2A'^^)] and G = VAB/^, for A = -trE^ and 7 = ^Eg : Eg. If trE'' > 0, 
both K and G are 0. 



Plastic flow rule: 



where P -tr(M)/3, r = ^Mo : Mo/2, and ^ = t/P. If /i < then = 0. 



As for the initial conditions, we declare that the granular body begins free of any plastic defor- 
mation, and that the body's initial deformation increment is elastic: 



The granular model being proposed in this paper is fully defined by the closed system of equations 
H [m [H [H [m [H [13 and [II] along with these initial conditions. 

6. Implementation and results 

The above equations were solved numerically using the ABAQUS/Explicit software package 
(from ABAQUS 6.5) with the constitutive behavior coded as a Vectorized User Material (VUMAT) 
model. In line with our interests, we choose three standard geometries that eventually have a 
sustained, well-developed flow response: 

• Long inclined chute 

• Annular Couette cell (with downward gravity) 

• Wide flat-bottom silo flow 

The simulations are run long enough for transient behavior to vanish. We may refer to this 
as "steady" flow, but a more general meaning is implied, since environments like silo drainage do 
not correspond to an Eulerian boundary value problem (the top surface always descends). Beyond 
visual observation, the disappearance of transients is quantiflable in the total kinetic energy, which 
becomes constant (compared to the kinetic energy transients) when flow is well-developed. 

6.1. Numerical considerations: The method of inertial density reduction 

ABAQUS 6.5 implements Lagrangian flnite element deformation, presenting some challenges in 
the modeling of steady flow behavior. In particular, there are two major numerical concerns that 
shall now be addressed: 

1. The steady response generally emerges after a high magnitude of deformation. However, accu- 
racy decreases dramatically when element distortion is too large. Arbitrary Lagrange-Eulerian 
(ALE) routines can be used to counteract this by periodically sweeping the mesh to a less dis- 
torted configuration. However, we find that the convection of variables between mesh sweeps 
carries unacceptably large error especially near boundaries, where spatial gradients have a re- 
duced sample space. Later versions of ABAQUS define purely Eulerian elements, which may 
improve this situation. However, as will be discussed momentarily, a different remedy can be 
instituted that uses the typical formulation. 




(18) 



FP{t = 0) = 1, LP(t = 0) = 0. 
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2. The constitutive relations are increasingly sensitive in the low pressure regime. By equation 
[T8l low pressure material at the onset of yield plastically flows at a rate proportional to P~^/^, 
which can cause rapid variation in during the transient phase. From an implementation 
standpoint this induces numerical stiffness; the rate of change may be too large to be adequately 
represented by the stable time increment. Furthermore, by equation 1171 the elastic moduli 
vanish sharply with decreasing pressure: k oc P^/^ as P — s- 0. Once again, this requires that 
the time step be small enough at low pressures so as to accurately track the changes in elastic 
properties. Numerical stiffness of the type just described is an important consideration since 
several of the flow environments to be investigated utilize the traction-free boundary condition, 
inducing non-negligible regions of low pressure in the material body. 

A single encomassing remedy, novel to this author's knowledge, was instituted to resolve both 
these issues: the inertial density of the material was artificially reduced by several orders of mag- 
nitude, holding the gravitational density fixed. The true granular material density (roughly 1500 
kg/m'^ for glass beads) was used in body force computations in order to accurately represent the 
force of gravity. However, the density used in computing the inertial force pv was artificially de- 
creased by factors on the order of 10^ — 10^. In systems without excessive curvature, the inertial 
force should vanish in steady state, so this alteration does not impact the correctness of the well- 
developed flow solution. To be clear, density reduction is performed only on the unstressed reference 
body, so that mass conservation is upheld during flow. 

By decreasing the inertial density, the wave speed increases, causing transients to pass quickly 
in physical time. This allows the steady state velocity and stress fields to emerge faster than 
the time necessary for large distortions to occur to the body (and likewise the mesh). A clear 
analogy for this effect can be observed in a rudimentary mass/spring/damper system in series — 
under applied force, the distance over which the motion of the mass is unsteady decreases as mass 
decreases. Using this technique, most simulations of the granular model reached steady state before 
any nodal displacements were noticeable to the eye. At the same time, since the stable time 
increment is proportional to y^, high enough density reduction reduces the time step (relative to 
the deformation speed) below the threshold necessary to accurately represent low pressure material 
during the transient phase. Hence, artificial density reduction confronts both issues raised above. It 
should be noted though, that there is no apparent benefit in terms of computation time when using 
this method since the transient time period and the stable time increment both decrease. More 
details on the method can be found in [46| . 

6.2. Verifying the numerical model 

Throughout this work, the model's six parameters are assigned the following values: 



B = 7 X 10'^ Pa 


lo = 0.279 


Ps = 2450 kg/m^ 


d = 0.003 m 


fis = tan(20.9°) 


fi2 = tan(32.76°) 



Recall that ^ is also a parameter, but its value is tied directly to the value /i2. Except for d, 
these values were all lifted from Jiang and Liu [i^l and Jop et al. As both groups considered 
spherical glass beads, it is assumed that their data should be representative of the same material. 
The particle diameter is set to 3mm, as this is the common value used in the experiments and DEM 
simulations of the MIT Dry Fluids Group, whose data this model will be compared against later. 
It should be noted that the particle diameter easily scales out in the non-dimensionalization. The 
elements used, unless otherwise stated, are hexahedral of type C3D8R in ABAQUS/Explicit. 

To check the implementation of the material model, one element tests were performed wherein 
a box element is compressed laterally with a fixed pressure, and sheared tangentially with various 
shear tractions. Implicit and explicit numerical integration routines were encoded into VUMAT 
files, and it was verified that both give steady fiow rates that match the analytical form (equation 
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[7|) to high accuracy. Since exphcit integration is only conditionally stable, the time step must be 
sufficiently small in order for numerical stability at the higher shear stress values. The implicit 
routine is always stable numerically, but for higher shear stresses, the time step must be reduced to 
ensure convergence of the Newton-Raphson solver. For the simulations detailed in this work, the 
time step is always small enough for explicit stability, a byproduct of the artificial density reduction. 
Together with the fact that explicit integration is computationally faster, the choice was made in 
the upcoming simulations to use the explicit routine. 

6.3. Rough-walled inclined chute 

Figure [H^a) reviews the geometry and boundary conditions for the rough inclined chute envi- 
ronment, which was originally studied by Jop et ai. In that work, equation [5] was directly solved 
as a non-Newtonian Bingham fluid, letting viscosity go to oo below yield. The approach delivers a 
steady velocity computation (with some numerical difficulty near the rigid interface), but by avoid- 
ing elasticity altogether, no stresses can be computed in the rigid zones. In terms of completeness, 
the Bingham fluid approach is useful, but not equipped to handle the general traction/kinematic 
boundary value problem. The elasto-plastic formulation herein is complete but has added many 
new relationships and a different kinematical perspective on deformation. However, the predicted 
steady flow should not differ noticeably from the Bingham flow, as elastic stretching vanishes in 
steady-state. By analyzing this case, we hope to show sufficient agreement between the two solu- 
tions, ensuring the elasto-plastic model has inherited the successes of the Bingham flow solution in 
this environment [i^ . 

The chute is long such that the flow can be deemed invariant in the z direction. The bottom of 
the chute is at y = 60d and the chute is 142c? wide. The angle of incline of the chute with respect 
to the horizontal is 22.6°. 

The chute model uses 2600 box elements arranged in a grid 65 elements wide in the x direction 
by 40 in the y. The assembly is only 1 element wide in the z direction — symmetry conditions 
were invoked by constraining the nodes on the back face of the xy slab to move with the same 
displacements as their front face counterparts. The simulation must start in a reference state with 
no gravity. Gravity is then ramped up to its full value. Since the material must be compressed 
before any shear stress can be supported, the components of gravity are applied one at a time; the 
y body force is first smoothly ramped up over 1 x 10~^s to its final value of cos(22.6°), where 
fg = 0.63gps. After a 2 x 10~^s pause for any needed relaxation, the component pointing down the 
chute is smoothly ramped up to fg sin(22.6°) over 2 x lO^^s. The assembly is then left to fiow until 
a total time of 5 x 10~'*s has elapsed. 

The inertial density was artificially decreased by a factor of several hundred thousand to a value 
of p = 2 X 10^'^. At the free surface, theoretically, the pressure goes to zero causing the elastic 
moduli to vanish. This is dangerous for procedures such as ABAQUS/Explicit that consistently 
send waves through the material — any small perturbation could cause a node to accelerate out of 
control. Moreover, the flow rule becomes undefined. To keep the free surface in tact, a few Pascals 
of overpressure are applied. 

While the flow is eventually very steady to the eye, a more quantitative measure is desired. 
From t = — 4 X 10~^s, the motion is markedly transient due to gravity ramping up. The system 
then relaxes toward the steady state, first rapidly but then slowly as steady state is approached. 
Comparing the system's rate of relative kinetic energy change in the fast relaxation period to its 
value at the end gives a rough criterion for how steady a flow is. At the end of the run, the rate of 
decrease of the relative kinetic energy is over 500 times smaller than the value it takes during the 
initial relaxation, which indicates a strong steady-state behavior. 

Comparisons to the numerical results of Jop et al. are displayed in flgurelUb). The agreement 
is quite good considering how disparate the solving methods are. Jop utilized a fixed grid finite- 
difference scheme solving a non-Newtonian Navier-Stokes-type problem, while the elasto-plastic 
results were obtained by a Lagrangian explicit procedure on nodes that are constantly moving 
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Fully rough walls Free surface 




Figure 4: (a) The rough inchned chute setup, (b) The velocity as a function of x at depths y/d = 0,9.2,18.4,27.6. 
The dashed curves are the numerical results of Jop et al. [44| via finite differences with a Stokcs-type solver. The 
solid curves are from the steady solution of the elasto-plastic model implemented on ABAQUS/Explicit. (c)-(d) 
Elasto-plastic results viewing the xy plane with y downward: In (c), regions of plastic flow arc in red and static in 
blue. In (d), the normalized pressure P/pgd is displayed. 

during steady-state. While differences of numerical procedure can cause different sources of gain 
and loss, it seems more likely that the discrepancies are stemming from the free surface condition. 
The true deformation rate at the top corners is actually infinite. As this is numerically impossible 
for either scheme, large but finite gradients occur there numerically as determined by the choice of 
the free surface treatment. This has a clear trickle-down effect on the global velocity field and, in 
lieu of the particular fashion in which the solutions differ, we suspect this effect causes much of the 
difference between the two numerical data sets. 

One distinguishing feature of the elasto-plastic model is its ability to calculate both flow and 
stress everywhere. Figure lljd) displays the pressure distribution over the full geometry. Note 
that the stresses vary smoothly through the transition from yielding to static, as opposed to the 
interfacial stress issues that commonly occur with Bingham fluid models. In particular, observe 
that the pressure field goes from decreasing linearly in the flowing zone, to decreasing somewhat 
nonlinearly in the static zone. In this geometry, compressive stresses in the x, y and z directions of 
the steady flowing zone must all be identical under codirectionality. A hydrostatic pressure profile 
is induced as a result. But upon descending below the flowing layer, codirectionality no longer has 
this influence on the stresses and a somewhat more complicated elasto-static form for the pressure 
field becomes apparent. Discrete simulation data for this environment would be helpful for checking 
the validity of the stress profile in the static zone. Notably, the primitive fully-rough boundary 
conditions used may not be the most accurate reflection of the true conditions on the bottom and 
side walls. 

6.4- Annular Couette cell 

We now move on to the annular Couette cell, an environment which has no previously known 
solution to the Jop et al. flow law. The results shall be compared directly against the myriad of data 
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on this environment compiled by G. D. R. Midi [59| and thus the geometric specifics and boundary 
conditions were chosen so as to give a good representation of the conditions used in these studies. 
Referring to figure [HJa) for general details on the environment, the following values were selected: 
Wwaii = 1.25rad/s « 0.2rev/s, the distance from inner to outer wall is Xout = 30(i, the height is 
Zbottom = lOrf, and the inner wall radius is 40(i. At the walls, the material motion must match the 
wall motion in the 9 and x directions, but material can slide without resistance up and down the 
walls. 

Since fiow and stress should be symmetric in the direction, the behavior as seen in a downward 
cut through the annular trough should represent the global behavior. A narrow sector of the annulus 
(total angle 0.1°) was likewise simulated using periodic boundary conditions on the front and back 
faces — nodal displacements on the front face are constrained to be identical to those on the back face 
except rotated appropriately by 0.1°. The sector is modeled using 40 elements in the x direction, 
15 in the z, and a thickness of one element in the 9 direction, for a grand total of 600 elements. 
Almost all the motion is known to occur near the inner wall in this environment. After preliminary 
tests of the elasto-plastic model produced the same conclusion, a bias was utilized to crowd nodes 
closer to the inner wall and improve precision. The bias resulted in half of the elements occurring 
within a distance of 6d from the inner wall. 

The inertial density was reduced by a factor of over one hundred thousand to a value of p = 0.01. 
First, gravity is smoothly ramped up to its full value over < = — 1 x lO^^s. During this time 
period, as before, a slight overpressure is also applied to the free surface for stability. The rotation 
of the inner wall then commences, smoothly ramping up to a final value of LHwaii = 1.25rad/s over 
t = 1 — 2.5 X 10~^s. The simulation is then left to flow until a total time oit — \ x 5^'*s has elapsed. 

Once the inner wall rotation reaches ujwaii, the flow soon after begins to relax toward a macro- 
scopically visible steady-state. Evidencing the fact that the system has indeed reached sufficient 
steady-state by the end of the run, the rate of kinetic energy loss of the system at the end of the 
simulation is over three orders smaller than the loss rate immediately after the ramp-up. 

Observing figures O^b) andO^c), we notice a few major qualitative points. For one, the fiow 
forms a clear shear band near the inner cylinder, a fact corroborated by the provided data. While 
the predicted thickness is somewhat smaller than found experimentally, the fact that it is on the 
same order is a major highlight for a model with no added fitting parameters (though we take this 
point with a grain of caution as shall be discussed momentarily). The authors of [44| expressed 
their belief that the flow rule would be incapable of describing narrow shear bands. However, the 
above results show quite clearly that this is not the case for the elasto-plastic model. 

Another observation is that the velocity profile does not vary to any observable extent in the z 
direction. This result has been verified in DEM simulations of this environment [48i] . where it was 
found that almost no vertical fluctuation occurs in the fast zone near the inner wall. Of course, 
the elasto-plastic velocity field dies off much differently than experiment. Where the experimental 
data is shown to be well-fit with an exponential decay that extends throughout the Couette cell, the 
elasto-plastic solution predicts a sharp cutoff around 5.5cJ from the inner wall. This result should 
come as no surprise, since the model does not account for quasi-static, non-local behavior, of which 
slow exponential decay is a textbook case. 

Viewing figure [S]Jd), it is clear that while the velocity appears invariant in the z direction, other 
important stress-based quantities are not. The stress ratio [i = t jP indeed has a highly non-trivial 
spatial distribution, refiective of the nonlinear dependence of the stress state on the plastic flow. 

Experimental data on this flow environment indicate that the shear band width stays relatively 
constant for a range of slower wall speeds |59i] , and grows with increasing wall velocity in the faster 
flow regime ^49i] . As expected from rate-sensitivity, we have verified with multiple simulations that 
the elasto-plastic model agrees with the latter case: the shear band width correlates positively with 
inner wall speed (holding gravity constant) . We expect this correlation to falsely persist even as the 
wall speed becomes small. This point serves as a careful reminder that the model has no internal 
length-scales and requires user discretion when assessing flow predictions — one must check that the 
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Figure 5: (a) The annular Couette setup, (b) The velocity profile normalized by the wall speed as predicted by 
the elasto-plastic model at half-height (thick green line) compared against fifteen experimental and discrete element 
simulation data sets for this type of fiow as compiled in 59]. In (c)-{d) the xz-plane is shown with z downward: (c) 
The velocity normalized by the wall speed, (d) The stress ratio /x. 

normalized flow rate and shear-rate gradients obey the limitations outlined in sections 14.21 and [2l 



6.5. Flat-bottomed silo 

While the past two environments have focused on comparing velocity profiles, the stresses have 
yet to be directly tested. Experimentally, the stress tensor is a difficult quantity to measure. The 
stresses in 2D disk assemblies can be approximated using photoelastic grain material, however, there 
is not currently an experimental method available to measure the stress tensor within an arbitrarily 
flowing 3D granular material. For this measurement, the best option as yet is to utilize DEM and 
compute the local stress tensor from grain fluctuations and inter-grain contact forces per equation 

M 

Rycroft et al. [t^ have performed DEM simulations of wide silo drainage. The computed stress 
and flow fields shall now be compared directly against the predictions of the elasto-plastic model. 
A schematic diagram of the flow geometry is pictured in figure [Gja). In accordance with Rycroft's 
simulations, we model the silo as having an opening width of Qd, a height of TOd, and a full width of 
150(i. Rycroft enforced periodic boundary conditions in the z direction giving the silo an apparent 
z thickness of but without wall-ordering effects. Since the flow should not fluctuate in the z 
direction, we simply enforce plane-strain conditions (using plane-strain elements of type CPE4R). 
Furthermore, the silo has left-right symmetry about the vertical center-line, which we take advantage 
of by modeling only the right half of the silo. 

The floor of the silo is modeled as having a frictional interaction with the material character- 
ized by a coefficient of friction /ifioor = 0.2. This number was estimated loosely from Rycroft's 
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Figure 6: (a) The flat-bottomed silo setup. In (b)-(c), a comparison of elasto-plastic velocity results ( — ) to DEM data 
(- -) c/o Rycroft [76[ . (b) The y velocity component as a function of x at heights y = 5d, lOd, 30d. (c) Trajectories 
predicted by the elasto-plastic model alongside the DEM trajectories. (Container outline provided for ease of viewing.) 



simulation — to account for the interaction induced by the effects of particle rolling, rearranging, 
and dragging along the floor, the element/floor interaction was modeled by reducing Rycroft's 
floor/particle friction coefhcient by 60%. Future study would be necessary to determine the com- 
plete and precise form for element/floor sliding interactions. While Rycroft's simulation utilized side 
walls made of the same frictional material as the floor, for simplicity, the elasto-plastic simulation 
employs the simpler condition of no x displacement at the side-walls. This could be enhanced in 
the future, but in wide silos, the details of the side-walls have only a small effect on the dominating 
behavior. 

Ideally speaking, the boundary condition at the silo orifice should be zero stress tractions. 
However, this is highly problematic numerically. In reality, silo flows develop a "free-fall arch" 
directly above the orifice [s^l- The arch is typically only a few particles high and connects the 
edges of the opening. Once a particle passes through this hypothetical arch, it enters free-fall and 
becomes gas-like. A granular material element within the free-fall arch would realistically have a 
smaller packing fraction, but still support compressive stresses through internal random particle 
collisions. The elasto-plastic model, however, does not include gas-like effects. Such dilation would 
be assigned to the elastic deformation, causing the elastic moduli to vanish. 

Our interest is not in the details within the immediate vicinity of the orifice, rather the bulk 
material behavior within the greater silo apparatus. However, with a zero-traction condition at the 
opening, the situation described above destabilizes the simulation prematurely. While free surfaces 
are usually taken care of by adding a slight amount of artificial compression, this remedy will not 
suffice here because the applied pressure has too much of an effect on the evolution of the outflow 
rate. 



19 



We are left with the alternative of using kinematic boundary conditions at the orifice. It would 
be overreaching to assign any particular velocity profile at the orifice. Instead, we fix the total flux 
out the orifice and let the material response determine the shape of the flow profile. To match the 
outfl-ux in Rycroft's simulation, 



was instituted at the orifice, encoded as an equation constraint in ABAQUS. 

Not far from the opening, large inhomogeneous deformation occurs at small length-scales, ne- 
cessitating many small elements to maintain accuracy. A grand total of 9750 elements were used in 
modeling the half-silo. To minimize discretization error, the orifice was modeled with a half-width 
of 15 elements. The adjacent silo floor was modeled 60 elements wide. The silo height was modeled 
with 130 elements. The element width was constant within the orifice, but bias was used along 
the other boundaries to maintain smooth changes in element sizes throughout. Elements shrink 
vertically as a sole function of distance from the silo bottom. The element width is uniform at 
the top surface, but moving downward, elements crowd the center so the floor region has a smooth 
transition from wider elements near the wall to narrower elements adjacent to the orifice. 

Due to the high number of elements and small minimal element size to system size ratio, this 
flow is too computationally intensive for one one processor. Instead, the 12- node Truesdell cluster of 
the MIT Solid Mechanics Group was employed to solve the problem in parallel. Using domain-level 
parallelization, the cluster split the half-silo into 12 spatial regions. 

The problem solved was as follows: From i = — 5 x 10~^s, gravity is gradually turned on 
while constraining the nodes along the opening from moving in the y direction in order to model a 
closed orifice. Then, from t = 0.55 — 5.00 x 10~^s, the orifice is gradually opened by ramping up the 
outflow rate from zero to Qout/2. The simulation is then left to flow until a total time oit — 10^'^s 
has been reached. 

This environment does not possess a steady-state since it is not an Eulerian boundary-value 
problem at the free surface. Instead, patternistic behavior eventually occurs, which signifies that 
transients have finished passing — starting at approximately t ~ lO^^s, the velocity and dP fields 
appear to fluctuate regularly. 

The first direct comparison that should be made is between the elasto-plastic and DEM flow 
profiles. To represent fully-developed mean behavior, Rycroft's flow data was averaged over 100 
frames, during a period of what appears to be transient-free flow. Similarly, the elasto-plastic flow 
was also averaged over many frames of fully-developed motion. To improve the validity of this 
average, the model was run an extra 5 x 10~'*s longer and the time average was performed over the 
range t = 0.5 — 1.5 x 10~'*s, comprising 127 frames. Figure [6] displays the comparison. Overall the 
agreement is sufficient. The particular way in which the peak in the downward velocity component 
broadens as height increases is well captured by the model. Once again, as is now a common theme, 
the elasto-plastic downward velocity appears to change more rapidly in space. Non-local effects 
such as diffusion could smooth out these sharper variations and possibly improve the agreement. 
Observing figure [G^c), the DEM and elasto-plastic trajectories agree well, especially below y k, A{)d. 
The differences in the upper silo trajectories could stem from the fact that in the DEM simulation, 
the top free surface lowered notably by the time steady flow could be imaged, whereas the elasto- 
plastic simulation reached this point before any noticeable drop of the top free surface (due to the 
artificial density reduction). 

To compare instantaneous behavior, we observe snapshots of the shearing rate profile. As shown 
in figure [71 well-developed elasto-plastic behavior involves two long "arms" of shearing that extend 
from the edges of the orifice to the top surface. From those arms, narrower, weaker shear bands 
(note the log-scale) continually drop down, creating an intersecting "mesh". The DEM solution 
does look similar, with two dominant arms of shearing. Note as well that near the top of the silo, 
both plots show the shearing spreads out approaching the free surface. The most obvious difference. 




Right half-opening 



Vy{x,y^ 0)dx = 2.19 x 10"^ — 
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Figure 7: The plastic shearing rate dP, expressed in units of gjd, plotted in logarithm form to accentuate small 
features, (a) The elasto-plastic solution: Note the intricate pattern of shear bands that fill the region between the 
two long shearing arms. The long-time behavior has the bands fall down one-by-one from the larger shearing arms, 
(b) The DEM solution c/o Rycroft [T^I : Similar to the elasto-plastic except blurred out by the box-averaging and the 
effects of non-local "diffusion" that are ignored by first-order elasto-plasticity. Both plots use the same color scale. 
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Figure 8: The instantaneous deviatoric principal stress directions plotted as crosshairs with lengths corresponding to 
the associated deviatoric stress eigenvalues, (a) The elasto-plastic solution: The thicker of the lines corresponds to 
the major principal stress direction, (b) The DEM solution c/o Rycroft l76il : Major principal stress in purple, minor 
(and intermediate where visible) in blue. 
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as has been previously noted, is that the elasto-plastic solution has sharp flow cutoffs whereas the 
shear rate always gradually tails off in the DEM. To reiterate, the model lacks non-local quasi-static 
terms so flow cutoffs are to be expected. Even if narrow bands did exist in the DEM, the box-average 
being performed would most assuredly obscure them. 

Moving on to the stresses, we first check the principal stresses and directions during fully- 
developed flow. As is evident in figure [51 the principal stress orientations predicted by the elasto- 
plastic model closely match those of the DEM. Both show the major principal stresses forming 
"arches" about the orifice and, moving away from the opening, the principal stress chains transition 
to becoming more vertical. In the DEM, it is apparent that the major principal stresses adjacent to 
the side walls have a slightly tilted orientation, whereas those of the elasto-plastic solution remain 
almost perfectly vertical. This is entirely due to the chosen boundary conditions. The DEM utilized 
frictional sidewalls and the tilt indicates the walls are exerting an upward shear on the material. 
This wall shear can be traced back to the filling process, where pouring causes an "active wall state" 
that resists the downward motion of the grains. On the contrary, the side walls in the elasto-plastic 
simulation are frictionless. A side wall friction may be included in the future to better model this 
minor effect. 

For a more direct comparison of the relative sizes of the stress invariants, see figure (5] Qualita- 
tively, the two solutions show similar spatial changes. Both show a region of lower /j, that swoops up 
from the lower side walls together with scattered minima of fi in the upper-middle region of the silo. 
Quantitatively, it appears the fj. values in the flowing zone are lower for the elasto-plastic solution. 
A probable explanation for this could be that the grains in the DEM simulation have a different 
surface roughness than the glass beads of [i^l and should in actuality have higher /is and fj,2 values. 
Recall that our material parameters were not extracted from the DEM data, but come from two 
different papers. In [l7| it was demonstrated that increasing the surface roughness of grains causes 
Us to increase. If the elasto-plastic model were implemented with a larger friction constants, the ^ 
proflle would indeed increase in the flowing regions, as larger fi would be needed to invoke the same 
plastic shear rate. 



7. Conclusion and future directions 

This work has demonstrated a highly general, 3D granular continuum model that unifles the re- 
cent results in granular elasticity and plasticity. The unification follows a rigorous F'^F^ decomposi- 
tion, with elastic and plastic response combined using flnitc-dcformation elasto-plasticity. The model 
can be used to predict flows uniquely in any environment with mechanically well-posed boundary 
conditions and/or body forces. The model was implemented as a user-material in AB AQUS/Explicit 
and tested in three unrelated geometries. With no fltting, it appears to give qualitative, and in some 
case quantitative, predictions for both the stress and velocity fleld in arbitrary granular flow geome- 
tries. Even so, there are a few clear avenues of future work with regard to improving the current 
model. 

7.1. Quasi-static non-locality 

The most glaring effect absent from the model is that it cannot account for "blurring" in the 
flow flelds. As described in depth in sections |4. 21 and when the normalized shear rate decreases 
low enough, the rheology is no longer determinable from a simple relation of the form / = g~^(fi). 

One would speculate a more complete form has higher order gradients in stress, flow rate, and/or 
state parameters balanced by some additional grain-level length-scale. In regions of moderate flow 
rate, these spatially second order effects should be dominated by the local rheology. This is clear 
from the above tests, which show the model does indeed perform better where the flow is faster. 
But equally clear is the fact that these terms cannot be ignored near static yield. Section 14.31 
lists several possible theories to describe non-local behavior, but how each candidate would fit 
theoretically within the current model remains to be studied. 
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By neglecting quasi-static non-locality, the model is tacitly assuming that the material is ca- 
pable of forming clean solid-like/fluid-like interfaces. This is almost always an idealization that 
is unrealized, as particles in a random packing rarely have the geometric ability to assemble in a 
fashion that avoids overlap with the predicted interface. A particle on the interface, being unable to 
both shear with the fluid and remain static with the solid, instead transmits some of the shearing 
behavior from the fluid-like zone into the solid zone, thereby explaining the gradual tails we see in 
actual granular flow profiles. Under particular circumstances, however, granular flow can be made 
to segregate clearly into flowing and completely static zones. Thompson and Grest [s^ have shown 
that a monodisperse 2D disk assembly undergoing horizontal planar shear with downward gravity 
does indeed have a flow cutoff with zero flow occurring beneath a shear band at the top. This 
behavior is a rarity brought on by the fact that solid-like material in this geometry can and does 
form a hexagonal crystal and the horizontal fluid/solid interface happens to align perfectly with a 
crystal plane. 

1.2. Dilation 

The present model avoids all inelastic dilation. Though an argument for why this is acceptable 
for our current purposes is presented in section 14. 4| certain benefits would come with properly 
accounting for the small amount of dilation that occurs in dense flow. 

Inelastic dilation comes in two primary forms: plastic and gas-like. Shear and shearing dilation, 
as described in section 14.21 are prototypical plastic mechanisms. Gas-like dilation can occur in 
particularly energetic surroundings, even below the plastic yield criterion — for example, consider 
shaking a box of rubber spheres under uniform pressure. Granted, gas-like behavior is outside the 
dense flow regime we currently study. Even so, it is possible for a small region of an otherwise dense 
flow to become gas-like, and consequently a means of dealing with this behavior is desired. Recall 
this was a crucial issue in the silo geometry, which is almost entirely dense except within the small 
free-fall arch that encompasses the opening. Gas-like effects may require a granular temperature 
and heat flux to express the internal pressure. 

Several models for plastic dilation have been proposed (see section l44l) . but a direct 3D discrete 
element study would be ideal for quantifying the precise dilatational dependences. One might 
hypothesize, based on equation [Sj that the plastic dilation has a form 

^ = -J^ = A{t^)B{^{I) - 0.636-") 
d-fP dP{M) ^" ' ' 

where t] = log(detFP) measures plastic dilation, 7^ is a plastic shear strain, / = tF(M)/ {^\J^^ is 
the inertial number, and the functions A and B are empirical, with B = 1 when the magnitude of its 
argument is large and B ^ as the input goes to zero. The above states that a material originally 
at random close packing ($ = 0.63) obeys Bagnold type dependences at steady flow. While the 
flow is unsteady, it dilates according to A and then ceases dilation as the Bagnold relationship 
is approached. A relation of this form appears to agree with results of Rycroft et al. [76], and 
collaborative efforts with Rycroft are underway to quantify this evolution law. 

Once a form for the plastic dilation has been verified, it would serve to enhance the computation 
of elasto-statics, where the moduli are known to vary with the packing fraction, as well as give 
meaningful packing fraction data throughout. However, the technique of artificial density reduction 
may no longer be valid, since the steady packing fraction, while only a few percent different than 
at the start, may take more time to develop than is typically allotted in a simulation run. 

7.3. Flow condition 

The plastic flow rule of Jop et al. asserts the codirectionality condition. A direct test of the 
flow condition should be performed to verify whether this is indeed the best candidate. Such a 
test may be forthcoming with recent DEM data on large conical hopper flow. Some models have 
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had success utilizing double-shearing conditions and Rycroft's silo flow data has shown some 
deviation from codirectional flow. Even so, codirectionality has fit the current needs and appears to 
be sufficient to predict basic flow behavior. A fine-tuned flow condition would be essential to model 
highly asymmetric 3D flows. 
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